home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
PD Collection CD 1
/
PD Collection CD 1.iso
/
programer2
/
pari2
/
pari
/
c
/
alglin2
< prev
next >
Wrap
Text File
|
1991-12-08
|
34KB
|
1,362 lines
/********************************************************************/
/********************************************************************/
/** **/
/** ++++++++++++++++++++++++++++++ **/
/** + + **/
/** + ALGEBRE LINEAIRE + **/
/** + + **/
/** ++++++++++++++++++++++++++++++ **/
/** **/
/** (deuxieme partie) **/
/** **/
/** copyright Babe Cool **/
/** **/
/** **/
/********************************************************************/
/********************************************************************/
# include "genpari.h"
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* POLYNOME CARACTERISTIQUE */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN caract(x,v)
GEN x;
int v;
{
long n,k,l,tetpil,tx=typ(x);
GEN p1,p2,p3,p4,p5,p6;
switch(tx)
{
case 1: case 2: case 3: case 4: case 5: case 7: p1=cgetg(4,10);
p1[1]=0x1000004+(v<<16);p1[2]=lneg(x);p1[3]=un;return p1;
case 6: case 8: p1=cgetg(5,10);p1[1]=0x1000005+(v<<16);p1[2]=lnorm(x);
l=avma;p2=trace(x[1]);tetpil=avma;p1[3]=lpile(l,tetpil,gneg(p2));
p1[4]=un;return p1;
case 9: l=avma;p1=gsub(polx[255],x[2]);tetpil=avma;
p1=subres(x[1],p1);
if(varn(p1)==255) {setvarn(p1,v);return gerepile(l,tetpil,p1);}
else {tetpil=avma;return gerepile(l,tetpil,gsubst(p1,255,polx[v]));}
case 19: n=lg(x)-1;if(!n) return polun[v];
if((n+1)!=lg(x[1])) err(mattype1);
l=avma;p1=gzero;p2=gun;
if(n%2) p2=gneg(p2);p5=cgetg(4,10);
p5[1]=0x01000004;p5[3]=un;setvarn(p5,v);
p6=cgeti(3);p5[2]=zero;
p6[1]=0xff000003;p4=cgetg(3,14);p4[2]=(long)p5;
for(k=0;k<=n;k++)
{
p3=det(gsub(gscalsmat(k,n),x));p4[1]=lmul(p3,p2);p6[2]=k;
p1=gadd(p4,p1);p5[2]=(long)p6;
if(k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);
}
p2=mpfact(n);tetpil=avma;return gerepile(l,tetpil,gdiv(p1[1],p2));
default: err(mattype1);
}
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* Methode des traces :
ce programme retourne le polynome caracteristique,
et si un pointeur non nul est fourni,celui pointe
sur la matrice adjointe a la sortie. */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN caradj(x,v,py)
GEN x,*py;
long v;
{
long i,j,k,l,t,av1,av2,av3,av4,decal;
GEN p,y,z;
if(typ(x)!=19) err(mattype1);
l=lg(x);
p=cgetg(l+2,10);setvarn(p,v);
p[l+1]=un;
av1=avma;t=ltrace(x);av2=avma;
t=lpile(av1,av2,gneg(t));
p[l]=t;
av1=avma;
y=cgetg(l,19);
for (i=1;i<l;i++) y[i]=lgetg(l,18);
for (i=1;i<l;i++)
for (j=1;j<l;j++)
{
if (i==j) coeff(y,i,j)=ladd(coeff(x,i,j),t);
else coeff(y,i,j)=coeff(x,i,j);
}
for (k=2;k<l-1;k++)
{
z=gmul(x,y);
t=ltrace(z);av2=avma;
t=ldivgs(t,-k);av3=avma;
y=cgetg(l,19);
for (i=1;i<l;i++) y[i]=lgetg(l,18);
for (i=1;i<l;i++) for (j=1;j<l;j++)
if (i==j) coeff(y,i,j)=ladd(coeff(z,i,j),t);
else coeff(y,i,j)=lcopy(coeff(z,i,j));
av4=avma;decal=lpile(av1,av2,0);
p[l-k+1]=adecaler(t,av2,av4)?t+decal:t;
if(adecaler(y,av2,av4)) y+=(decal>>2);
av1=av3+decal;
}
t=zero;
for (i=1;i<l;i++)
t=ladd(t,gmul(coeff(x,1,i),coeff(y,i,1)));
av2=avma;t=lneg(t);
if ((long) py)
{
z=cgetg(l,19);
for (i=1;i<l;i++) z[i]=lgetg(l,18);
for (i=1;i<l;i++) for (j=1;j<l;j++)
coeff(z,i,j)=lcopy(coeff(y,i,j));
av4=avma;decal=lpile(av1,av2,0);
p[2]=adecaler(t,av2,av4)?t+decal:t;
*py=adecaler(z,av2,av4)?z+(decal>>2):z;
}
else p[2]=lpile(av1,av2,t);
p[1]=0x01000000+l+2;
return p;
}
GEN adj(x)
GEN x;
{
GEN y;
caradj(x,255,&y);
return y;
}
GEN caradj0(x,v)
GEN x;
long v;
{
long tx=typ(x),l,tetpil;
GEN p1,p2;
switch(tx)
{
case 1: case 2: case 3: case 4: case 5: case 7: p1=cgetg(4,10);
p1[1]=0x1000004+(v<<16);p1[2]=lneg(x);p1[3]=un;return p1;
case 6: case 8: p1=cgetg(5,10);p1[1]=0x1000005+(v<<16);p1[2]=lnorm(x);
l=avma;p2=trace(x);tetpil=avma;p1[3]=lpile(l,tetpil,gneg(p2));
p1[4]=un;return p1;
case 9: l=avma;p1=gsub(polx[255],x[2]);tetpil=avma;
p1=subres(x[1],p1);
if(varn(p1)==255) {setvarn(p1,v);return gerepile(l,tetpil,p1);}
else {tetpil=avma;return gerepile(l,tetpil,gsubst(p1,255,polx[v]));}
default: return caradj(x,v,0);
}
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* INVERSION D'UNE MATRICE */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN invmulmat(a,b) /* calcule a^(-1)*b, b etant une matrice.
( Il faut : nblig(a)=nbcol(a)=nblig(b) ) */
GEN a,b;
{
long p,u,m,nbli,nbco,i,j,k,av,av1,av2,av3,av4;
GEN aa,x;
nbco=lg(b)-1;if(!nbco) return cgetg(1,19);
nbli=lg(a[1])-1;
if (lg(a)-1 != nbli) err(invmuler1);
if (nbli!=lg(b[1])-1) err(invmuler1);
av=avma;
x=cgetg(nbco+1,19);
for (j=1;j<=nbco;j++)
{
x[j]=lgetg(nbli+1,18);
for (i=1;i<=nbli;i++)
coeff(x,i,j)=coeff(b,i,j);
}
av1=avma;
aa=cgetg(nbli+1,19);
for (j=1;j<=nbli;j++)
{
aa[j]=lgetg(nbli+1,18);
for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
}
for (i=1;i<nbli;i++)
{
p=coeff(aa,i,i);k=i;
if (gcmp0(p))
{
for (k=i+1;(k<=nbli)&&gcmp0(coeff(aa,k,i));k++);
if (k>nbli) err(matinv2);
else
{
for (j=i;j<=nbli;j++)
{
u=coeff(aa,i,j);coeff(aa,i,j)=coeff(aa,k,j);
coeff(aa,k,j)=u;
}
for (j=1;j<=nbco;j++)
{
u=coeff(x,i,j);coeff(x,i,j)=coeff(x,k,j);
coeff(x,k,j)=u;
}
p=coeff(aa,i,i);
}
}
for (k=i+1;k<=nbli;k++)
{
m=coeff(aa,k,i);
if (!gcmp0(m))
{
m=ldiv(m,p);
for (j=i+1;j<=nbli;j++)
coeff(aa,k,j)=lsub(coeff(aa,k,j),gmul(m,coeff(aa,i,j)));
for (j=1;j<=nbco;j++)
coeff(x,k,j)=lsub(coeff(x,k,j),gmul(m,coeff(x,i,j)));
}
}
}
av2=avma;
p=coeff(aa,nbli,nbli);
if (gcmp0(p)) err(matinv2);
else
{
for (j=1;j<=nbco;j++)
{
coeff(x,nbli,j)=ldiv(coeff(x,nbli,j),p);
for (i=nbli-1;i>0;i--)
{
av3=avma;
m=coeff(x,i,j);
for (k=i+1;k<=nbli;k++)
m= lsub(m,gmul(coeff(aa,i,k),coeff(x,k,j)));
av4=avma;
coeff(x,i,j)=lpile(av3,av4,gdiv(m,coeff(aa,i,i)));
}
}
av=lpile(av1,av2,0);
for(i=1;i<=nbli;i++)
for(j=1;j<=nbco;j++)
if (coeff(x,i,j)<=av1) coeff(x,i,j)+=av;
}
return x;
}
GEN invmat(a)
GEN a;
{
long av=avma,tetpil;
GEN b;
b=gscalmat(un,lg(a)-1);tetpil=avma;
return gerepile(av,tetpil,invmulmat(a,b));
}
GEN invmulmatreel(a,b) /* calcule a^(-1)*b, b etant une matrice.
( Il faut : nblig(a)=nbcol(a)=nblig(b) ) */
GEN a,b;
{
long p,u,m,nbli,nbco,i,j,k,av,av1,av2,av3,av4;
GEN aa,x,p1;
nbco=lg(b)-1;nbli=lg(a[1])-1;
if (lg(a)-1 != nbli) err(invmuler1);
/* la verif nblig(b)=nblig(a) n'est pas faite ... */
av=avma;
x=cgetg(nbco+1,19);
for (j=1;j<=nbco;j++)
{
x[j]=lgetg(nbli+1,18);
for (i=1;i<=nbli;i++)
coeff(x,i,j)=coeff(b,i,j);
}
av1=avma;
aa=cgetg(nbli+1,19);
for (j=1;j<=nbli;j++)
{
aa[j]=lgetg(nbli+1,18);
for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
}
for (i=1;i<nbli;i++)
{
p=labs(coeff(aa,i,i));k=i;
for(j=i+1;j<=nbli;j++)
if(gcmp(p1=gabs(coeff(aa,j,i),3),p)>0) {p=(long)p1;k=j;}
if (gcmp0(p)) err(matinv2);
else
{
if(k>i)
{
for (j=i;j<=nbli;j++)
{
u=coeff(aa,i,j);coeff(aa,i,j)=coeff(aa,k,j);
coeff(aa,k,j)=u;
}
for (j=1;j<=nbco;j++)
{
u=coeff(x,i,j);coeff(x,i,j)=coeff(x,k,j);
coeff(x,k,j)=u;
}
}
p=coeff(aa,i,i);
}
for (k=i+1;k<=nbli;k++)
{
m=coeff(aa,k,i);
if (!gcmp0(m))
{
m=ldiv(m,p);
for (j=i+1;j<=nbli;j++)
coeff(aa,k,j)=lsub(coeff(aa,k,j),gmul(m,coeff(aa,i,j)));
for (j=1;j<=nbco;j++)
coeff(x,k,j)=lsub(coeff(x,k,j),gmul(m,coeff(x,i,j)));
}
}
}
av2=avma;
p=coeff(aa,nbli,nbli);
if (gcmp0(p)) err(matinv2);
else
{
for (j=1;j<=nbco;j++)
{
coeff(x,nbli,j)=ldiv(coeff(x,nbli,j),p);
for (i=nbli-1;i>0;i--)
{
av3=avma;
m=coeff(x,i,j);
for (k=i+1;k<=nbli;k++)
m= lsub(m,gmul(coeff(aa,i,k),coeff(x,k,j)));
av4=avma;
coeff(x,i,j)=lpile(av3,av4,gdiv(m,coeff(aa,i,i)));
}
}
av=lpile(av1,av2,0);
for(i=1;i<=nbli;i++)
for(j=1;j<=nbco;j++)
if (coeff(x,i,j)<=av1) coeff(x,i,j)+=av;
}
return x;
}
GEN invmatreel(a)
GEN a;
{
long av=avma,tetpil;
GEN b;
b=gscalmat(un,lg(a)-1);tetpil=avma;
return gerepile(av,tetpil,invmulmatreel(a,b));
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* FORME DE HESSENBERG D'UNE MATRICE */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN hess(x)
GEN x;
{
long tx=typ(x),lx=lg(x),av=avma,tetpil,m,i,j;
GEN p1,p2,p3,y;
if((tx!=19)||(lg(x[1])!=lx)) err(mattype1);
y=gcopy(x);
for(m=2;m<lx-1;m++)
{
p2=gzero;
for(i=m+1;(i<lx)&&(gcmp0(p2=(GEN)coeff(y,i,m-1)));i++);
if(!gcmp0(p2))
{
for(j=m-1;j<lx;j++)
{
p1=(GEN)coeff(y,i,j);coeff(y,i,j)=coeff(y,m,j);coeff(y,m,j)=(long)p1;
}
p1=(GEN)y[i];y[i]=y[m];y[m]=(long)p1;
for(i=m+1;i<lx;i++)
{
if(!gcmp0(p3=(GEN)coeff(y,i,m-1)))
{
p3=gdiv(p3,p2);coeff(y,i,m-1)=zero;
for(j=m;j<lx;j++)
coeff(y,i,j)=lsub(coeff(y,i,j),gmul(p3,coeff(y,m,j)));
for(j=1;j<lx;j++)
coeff(y,j,m)=ladd(coeff(y,j,m),gmul(p3,coeff(y,j,i)));
}
}
}
}
tetpil=avma;return gerepile(av,tetpil,gcopy(y));
}
GEN carhess(x,v)
GEN x;
long v;
{
long av=avma,tetpil,tx=typ(x),lx=lg(x),r,i;
GEN *y,p1,p2,p3,p4;
if((tx!=19)||(lg(x[1])!=lx)) err(mattype1);
y=(GEN *)newbloc(4*lx);
y[0]=polun[v];p1=hess(x);p2=polx[v];
for(r=1;r<lx;r++)
{
y[r]=gmul(y[r-1],gsub(p2,coeff(p1,r,r)));
p3=gun;p4=gzero;
for(i=1;i<r;i++)
{
p3=gmul(p3,coeff(p1,r-i+1,r-i));
p4=gadd(p4,gmul(gmul(p3,coeff(p1,r-i,r)),y[r-i-1]));
}
tetpil=avma;y[r]=gsub(y[r],p4);
}
p1=gerepile(av,tetpil,y[lx-1]);
killbloc(y);return p1;
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* NORME */
/* */
/* Cree un GEN pointant sur la norme de x */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN gnorm(x)
GEN x;
{
long l,tx,lx,i,tetpil;
GEN p1,p2,y;
switch(tx=typ(x))
{
case 1 :
case 2 : y=mpmul(x,x);break;
case 3 : err(normer1);
case 4 :
case 5 : y=gmul(x,x);break;
case 6 : l=avma;p1=gmul(x[1],x[1]);
p2=gmul(x[2],x[2]);
tetpil=avma;
y=gerepile(l,tetpil,gadd(p1,p2));break;
case 8 : l=avma;p1=(GEN)x[1];
if (gcmp0(p1[3]))
{
p2=gmul(p1[2],gmul(x[3],x[3]));
p1=gmul(x[2],x[2]);
tetpil=avma;
y=gerepile (l,tetpil,gadd(p1,p2));
}
else
{
p2=gmul(p1[2],gmul(x[3],x[3]));
p1=gmul(x[2],gadd(x[2],x[3]));
tetpil=avma;
y=gerepile(l ,tetpil,gadd(p1,p2));
}
break;
case 10:
case 11:
case 13:
case 14: l=avma;p1=gmul(gconj(x),x);tetpil=avma;
y=gerepile(l,tetpil,greal(p1));break;
case 9 : y=subres(x[1],x[2]);break;
case 17:
case 18:
case 19: lx=lg(x);y=cgetg(lx,tx);
for(i=1;i<lx;i++) y[i]=lnorm(x[i]);
break;
default: err(normer1);
}
return y;
}
GEN gnorml2(x)
GEN x;
{
GEN y,p1;
long i,tx=typ(x),lx=lg(x),av,tetpil;
if(tx<17) return gnorm(x);
y=gzero;
if(lx>1)
{
av=avma;y=gnorml2(x[1]);
for(i=2;i<lx;i++)
{
p1=gnorml2(x[i]);tetpil=avma;
y=gadd(p1,y);
}
if(lx>2) y=gerepile(av,tetpil,y);
}
return y;
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* CONJUGAISON */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN gconj(x)
GEN x;
{
long lx,i,tx=typ(x);
GEN z,p1;
switch(tx)
{
case 1 :
case 2 :
case 3 :
case 4 :
case 5 :
case 7 : z=gcopy(x);break;
case 6 : z=cgetg(3,6);z[1]=lcopy(x[1]);
z[2]=lneg(x[2]);
break;
case 8 : z=cgetg(4,8);z[1]=copyifstack(x[1]);
p1=(GEN)x[1];
z[3]=lneg(x[3]);
if(gcmp0(p1[3])) z[2]=lcopy(x[2]);
else z[2]=ladd(x[2],x[3]);
break;
case 10: lx=lg(x);z=cgetg(lx,tx);
z[1]=x[1];
for(i=2;i<lgef(x);i++)
z[i]=lconj(x[i]);
break;
case 11: lx=lg(x);z=cgetg(lx,tx);
z[1]=x[1];
if(!gcmp0(x))
{
for(i=2;i<lx;i++)
z[i]=lconj(x[i]);
}
break;
case 13:
case 14:
case 17:
case 18:
case 19: lx=lg(x);z=cgetg(lx,tx);
for(i=1;i<lx;i++)
z[i]=lconj(x[i]);
break;
default: err(conjer1);
}
return z;
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* PARTIES REELLE ET IMAGINAIRES */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN greal(x)
GEN x;
{
long lx,i,j,av,tetpil,tx=typ(x);
GEN p1,p2,z;
switch(tx)
{
case 1 :
case 2 :
case 4 :
case 5 : z=gcopy(x);break;
case 6 : z=gcopy(x[1]);break;
case 8 : z=gcopy(x[2]);break;
case 10: lx=lgef(x);av=avma;
for(i=lx-1;(i>=2)&&(gcmp0(greal(x[i])));i--);
avma=av;
if(i<2) {z=cgetg(2, tx);z[1]=2;}
else
{
z=cgetg(i+1,tx);z[1]=0x01000001+i;
for(j=2;j<=i;j++) z[j]=lreal(x[j]);
}
setvarn(z,varn(x));
break;
case 11: lx=lg(x);av=avma;
if(gcmp0(x)) {z=cgetg(2,tx);z[1]=x[1];}
else
{
for(i=2;(i<lx)&&(gcmp0(greal(x[i])));i++);
avma=av;
if(i==lx)
{
z=cgetg(2,tx); setvalp(z, lx - 2 + valp(x));
setsigne(z,0); setvarn(z, varn(x));
}
else
{
z=cgetg(lx-i+2,tx);
for(j = 2; j <= lx - i + 1; j++) z[j] = lreal(x[j + i - 2]);
z[1] = x[1]; setvalp(z, valp(x) + i - 2);
}
}
break;
case 13:
case 14: av=avma;p1=gadd(gmul(greal(x[1]),greal(x[2])),gmul(gimag(x[1]),gimag(x[2])));
p2=gadd(gsqr(greal(x[2])),gsqr(gimag(x[2])));tetpil=avma;
z=gerepile(av,tetpil,gdiv(p1,p2));break;
case 17:
case 18:
case 19: lx=lg(x);z=cgetg(lx,tx);
for(i=1;i<lx;i++) z[i]=lreal(x[i]);
break;
default: err(realer1);
}
return z;
}
GEN gimag(x)
GEN x;
{
long lx,i,j,av,tetpil,tx=typ(x);
GEN p1,p2,z;
switch(tx)
{
case 1 :
case 2 :
case 4 :
case 5 : z=gzero;break;
case 6 : z=gcopy(x[2]);
break;
case 8 : z=gcopy(x[3]);
break;
case 10: lx=lgef(x);av=avma;
for(i=lx-1;(i>=2)&&(gcmp0(gimag(x[i])));i--);
avma=av;
if(i<2) {z=cgetg(2, tx);z[1]=2;}
else
{
z=cgetg(i+1,tx);z[1]=0x01000001+i;
for(j=2;j<=i;j++) z[j]=limag(x[j]);
}
setvarn(z,varn(x));
break;
case 11: lx=lg(x);av=avma;
if(gcmp0(x)) {z=cgetg(2,tx);z[1]=x[1];}
else
{
for(i=2;(i<lx)&&(gcmp0(gimag(x[i])));i++);
avma=av;
if(i==lx)
{
z=cgetg(2,tx); setvalp(z, lx - 2 + valp(x));
setsigne(z,0); setvarn(z, varn(x));
}
else
{
z=cgetg(lx-i+2,tx);
for(j = 2; j <= lx - i + 1; j++) z[j] = limag(x[j + i - 2]);
z[1] = x[1]; setvalp(z, valp(x) + i - 2);
}
}
break;
case 13:
case 14: av=avma;p1=gsub(gmul(gimag(x[1]),greal(x[2])),gmul(greal(x[1]),gimag(x[2])));
p2=gadd(gsqr(greal(x[2])),gsqr(gimag(x[2])));tetpil=avma;
z=gerepile(av,tetpil,gdiv(p1,p2));break;
case 17:
case 18:
case 19: lx=lg(x);z=cgetg(lx,tx);
for(i=1;i<lx;i++)
z[i]=limag(x[i]);
break;
default: err(imager1);
}
return z;
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* */
/* TRACES */
/* */
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN assmat(x)
GEN x;
{
long lx,i,j;
GEN y,p1;
if((typ(x)!=10) || gcmp0(x)) err(mattype2);
lx=lgef(x)-2;y=cgetg(lx,19);
for(i=1;i<lx-1;i++)
{
p1=cgetg(lx,18);y[i]=(long)p1;
for(j=1;j<lx;j++)
{p1[j]=(j==(i+1)) ? un : zero;}
}
p1=cgetg(lx,18);y[i]=(long)p1;
if(gcmp1(x[lx+1]))
{
for(j=1;j<lx;j++)
p1[j]=lneg(x[j+1]);
}
else
{
gnegz(x[lx+1],x[lx+1]);
for(j=1;j<lx;j++)
p1[j]=ldiv(x[j+1],x[lx+1]);
gnegz(x[lx+1],x[lx+1]);
}
return y;
}
GEN trace(x)
GEN x;
{
long i,l,n,tx=typ(x),lx=lg(x),tetpil;
GEN y,p1,p2;
switch(tx)
{
case 1 :
case 2 :
case 4 :
case 5 : y=gmul2n(x,1);break;
case 6 : y=gmul2n(x[1],1);break;
case 8 : p1=(GEN)x[1];
if (!gcmp0(p1[3]))
{
l=avma;p2=gmul2n(x[2],1);
tetpil=avma;
y=gerepile(l,tetpil,gadd(x[3],p2));
}
else y=gmul2n(x[2],1);
break;
case 10: lx=lg(x);y=cgetg(lx,tx);
y[1]=x[1];
for(i=2;i<lgef(x);i++)
y[i]=ltrace(x[i]);
break;
case 11: lx=lg(x);y=cgetg(lx,tx);
y[1]=x[1];
if(!gcmp0(x))
{
for(i=2;i<lx;i++)
y[i]=ltrace(x[i]);
}
break;
case 9 : l=avma;p1=polsym(x[1],n=(lgef(x[1])-4));
p2=gzero;for(i=0;i<=n;i++) p2=gadd(p2,gmul(truecoeff(x[2],i),p1[i+1]));
tetpil=avma;y=gerepile(l,tetpil,gcopy(p2));
break;
case 13:
case 14: y=gadd(x,gconj(x));break;
case 17:
case 18: lx=lg(x);y=cgetg(lx,tx);
for(i=1;i<lx;i++)
y[i]=ltrace(x[i]);
break;
case 19: if (lx!=lg(x[1])) err(mattype3);
l=avma;p1=gcopy(coeff(x,1,1));
if(lx==2) return p1;
else
{
for(i=2;i<lx-1;i++)
p1=gadd(p1,coeff(x,i,i));
tetpil=avma;
y=gerepile(l,tetpil,gadd(p1,coeff(x,i,i)));
}
break;
default: err(mattype3);
}
return y;
}
/*===================================*/
/* Reduction en carres */
/*===================================*/
/*=======================================================
Reduction de Gauss ( Matrice definie >0 )
========================================================*/
GEN sqred1(a)
GEN a;
{
GEN b;
long av,av1,n,i,j,k,p,lim;
if (typ(a)!=19) err(kerer1);
if (lg(a[1])!=(n=lg(a))) err(mattype1);
lim=(avma+bot)>>1;
n--;av=avma;
b=gcopy(a);
for(i=1;i<=n;i++)
for(j=1;j<i;j++) coeff(b,i,j) = zero;
for(k=1;k<=n;k++)
{
if(gsigne(p=coeff(b,k,k))<=0) err(sqreder1);
for(i=k+1;i<=n;i++)
for(j=i;j<=n;j++)
coeff(b,i,j)=lsub(coeff(b,i,j),gdiv(gmul(coeff(b,k,i),coeff(b,k,j)),p));
for(j=k+1;j<=n;j++)
coeff(b,k,j)=ldiv(coeff(b,k,j),p);
if(avma<lim) {av1=avma;b=gerepile(av,av1,gcopy(b));}
}
av1=avma;
return gerepile(av,av1,gcopy(b));
}
GEN sqred3(a)
GEN a;
{
long n,av=avma,tetpil,lim,i,j,k,l;
GEN p1,z;
if (typ(a)!=19) err(kerer1);
if (lg(a[1])!=(n=lg(a))) err(mattype1);
lim=(avma+bot)>>1;
av=avma;z=cgetg(n,19);
for(j=1;j<n;j++)
{
p1=cgetg(n,18);z[j]=(long)p1;
for(i=1;i<n;i++) p1[i]=zero;
}
for(i=1;i<n;i++)
{
for(k=1;k<i;k++)
{
p1=gzero;
for(l=1;l<k;l++) p1=gadd(p1,gmul(gmul(coeff(z,l,l),coeff(z,k,l)),coeff(z,i,l)));
coeff(z,i,k)=ldiv(gsub(coeff(a,i,k),p1),coeff(z,k,k));
}
p1=gzero;
for(l=1;l<i;l++) p1=gadd(p1,gmul(gmul(coeff(z,l,l),coeff(z,i,l)),coeff(z,i,l)));
coeff(z,i,k)=lsub(coeff(a,i,i),p1);
if(avma<lim) {tetpil=avma;z=gerepile(av,tetpil,gcopy(z));}
}
tetpil=avma;return gerepile(av,tetpil,gcopy(z));
}
/*=======================================================
Reduction de Gauss (matrice symetrique quelconque)
Signature d'une matrice symetrique
( seule la partie superieure est consideree )
========================================================*/
GEN sqred2(a,flg)
GEN a; long flg;
{
GEN b,r;
long av,av1,av2,lim,n,i,j,k,l,p,sp,sn,t,u;
if (typ(a)!=19) err(kerer1);
if (lg(a[1])!=(n=lg(a))) err(mattype1);
av=avma;lim=(avma+bot)>>1;
r=cgeti(n);for(i=1;i<n;i++) r[i]=1;
av2=avma;b=gcopy(a);
t=(--n);sp=sn=0;
while (t)
{
for(k=1;(k<=n)&&(gcmp0(coeff(b,k,k))||(!r[k]));k++);
if(k<=n)
{
p=coeff(b,k,k);
if(signe(p)>0) sp++;else sn++;
r[k]=0;t--;
for(j=1;j<=n;j++)
coeff(b,k,j)=r[j] ? ldiv(coeff(b,k,j),p) : zero;
for(i=1;i<=n;i++) if (r[i])
{
for(j=1;j<=n;j++)
coeff(b,i,j)=r[j] ? lsub(coeff(b,i,j),gmul(gmul(coeff(b,k,i),coeff(b,k,j)),p)) : zero;
}
coeff(b,k,k)=p;
}
else
{
for(k=1;k<=n;k++) if (r[k])
{
for(l=k+1;(l<=n)&&(gcmp0(coeff(b,k,l))||(!r[l]));l++);
if(l<=n)
{
p=coeff(b,k,l);r[k]=r[l]=0;sp++;sn++;t-=2;
for(i=1;i<=n;i++) if(r[i])
{
for(j=1;j<=n;j++)
coeff(b,i,j)=r[j]? lsub(coeff(b,i,j),gdiv(gadd(gmul(coeff(b,k,i),coeff(b,l,j)),gmul(coeff(b,k,j),coeff(b,l,i))),p)) : zero;
}
for(j=1;j<=n;j++) if (r[j])
{
u=coeff(b,k,j);
coeff(b,k,j)=ldiv(gadd(u,coeff(b,l,j)),p);
coeff(b,l,j)=ldiv(gsub(u,coeff(b,l,j)),p);
}
coeff(b,k,l)=un;coeff(b,l,k)=lneg(un);
coeff(b,k,k)=lmul2n(p,-1);coeff(b,l,l)=lneg(coeff(b,k,k));
break;
}
if(avma<lim) {av1=avma;b=gerepile(av2,av1,gcopy(b));}
}
if(k>n) break;
}
}
if (flg) {av1=avma;return gerepile(av,av1,gcopy(b));}
else
{
avma=av;
b=cgetg(3,17);b[1]=lstoi(sp);b[2]=lstoi(sn);return b;
}
}
GEN sqred(a) { return sqred2(a,1); }
GEN signat(a) { return sqred2(a,0); }
/*===========================================================================
Diagonalisation d'une matrice symetrique REELLE;
matrice de passage orthogonale R
( Nouvelle version : 25/6/90 )
( Renvoie un vecteur a 2 comp :
1-re comp = vect des valeurs propres
2-me comp = matr des vecteurs propres ).
============================================================================*/
GEN jacobi(a,prec) GEN a;long prec;
{
long de,e,e1,e2,l,n,i,j,p,q,x,y,x1,y1,av1,av2,iter=0;
GEN c,s,t,u,ja,lambda,r,unr;
ja=cgetg(3,17);
n=(l=lg(a))-1;
e1=0x7fffff;
lambda=cgetg(l,18);ja[1]=(long)lambda;
for(j=1;j<=n;j++)
{
gaffect(coeff(a,j,j),x=lambda[j]=lgetr(prec));
if(((e=expo(x))<e1)&&(gsigne(x))) e1=e;
}
r=cgetg(l,19);ja[2]=(long)r;
for(j=1;j<=n;j++)
{
r[j]=lgetg(l,18);
for(i=1;i<l;i++)
affsr((i==j)? 1:0,coeff(r,i,j)=lgetr(prec));
}
av1=avma;
e2=expo(coeff(a,1,2));p=1;q=2;
c=cgetg(l,19);
for(j=1;j<=n;j++)
{
c[j]=lgetg(j,18);
for(i=1;i<j;i++)
{
gaffect(coeff(a,i,j),x=coeff(c,i,j)=lgetr(prec));
if((e=expo(x))>e2) {e2=e;p=i;q=j;}
}
}
a=c;
affsr(1,unr=cgetr(prec));
de=((prec-2)<<5);
while(e1-e2<de)
{
iter++;
/*calcul de la rotation associee dans le plan
des p et q-iemes vecteurs de base */
av2=avma;
x=ldivrr(subrr(lambda[q],lambda[p]),shiftr(coeff(a,p,q),1));
y=lmpsqrt(addrr(unr,mulrr(x,x)));
t=(gsigne(x)>0)? divrr(unr,addrr(x,y)) : divrr(unr,subrr(x,y));
c=divrr(unr,mpsqrt(addrr(unr,mulrr(t,t))));
s=mulrr(t,c);u=divrr(s,addrr(unr,c));
/* Recalcul des transformees successives de la matrice a et de la matrice
cumulee (r) des rotations : */
for(i=1;i<p;i++)
{
x=coeff(a,i,p); y=coeff(a,i,q);
x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
affrr(x1,coeff(a,i,p));affrr(y1,coeff(a,i,q));
}
for(i=p+1;i<q;i++)
{
x=coeff(a,p,i); y=coeff(a,i,q);
x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
affrr(x1,coeff(a,p,i));affrr(y1,coeff(a,i,q));
}
for(i=q+1;i<=n;i++)
{
x=coeff(a,p,i); y=coeff(a,q,i);
x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
affrr(x1,coeff(a,p,i));affrr(y1,coeff(a,q,i));
}
x=lambda[p]; y=coeff(a,p,q); subrrz(x,mulrr(t,y),lambda[p]);
x=y; y=lambda[q]; addrrz(y,mulrr(t,x),y);
/* if((e=expo(lambda[p]))<e1) e1=e;
if((e=expo(lambda[q]))<e1) e1=e; */
/* affsr(0,x); NON ! */
setexpo(x,expo(x)-de-1);
for(i=1;i<=n;i++)
{
x=coeff(r,i,p); y=coeff(r,i,q);
x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
affrr(x1,coeff(r,i,p));affrr(y1,coeff(r,i,q));
}
e2=expo(coeff(a,1,2));p=1;q=2;
for(j=1;j<=n;j++)
{
for(i=1;i<j;i++) if((e=expo(coeff(a,i,j)))>e2) {e2=e;p=i;q=j;}
for(i=j+1;i<=n;i++) if((e=expo(coeff(a,j,i)))>e2) {e2=e;p=j;q=i;}
}
avma=av2;
} /* Fin de la boucle (while) de recalcul */
avma=av1; return ja;
}
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~ ~*/
/*~ MATRICE RATIONNELLE-->ENTIERE ~*/
/*~ ~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
GEN matrixqz(x,pp)
GEN x,pp;
{
long av=avma,av1,tetpil,i,j,j1,m,n,t,fl,lim,nfact;
GEN p,p1,p2,p3,p4,p5,unmodp,pk,pt;
if(typ(x)!=19) err(matqzer1);
lim=(avma+bot)>>1;
n=lg(x)-1;if(!n) return gcopy(x);
m=lg(x[1])-1;if(n>m) err(matqzer2);
if(n==m)
{
p1=det(x);if(gcmp0(p1)) err(matqzer4);
return idmat(n);
}
p1=cgetg(n+1,19);
for(j=1;j<=n;j++)
{
p2=gun;p3=(GEN)x[j];
for(i=1;i<=m;i++)
{
t=typ(p3[i]);if((t>5)||(t==3)) err(matqzer3);
p2=ggcd(p2,p3[i]);
}
p1[j]=ldiv(p3,p2);
}
x=p1;
if(gcmp0(pp))
{
pt=gtrans(x);
p1=cgetg(n+1,19);
for(j=1;j<=n;j++) p1[j]=pt[j];
p2=det(p1);p1[n]=pt[n+1];p3=det(p1);
p4=ggcd(p2,p3);if(!signe(p4)) err(impl,"matrixqz when the first 2 dets are zero");
if(gcmp1(p4)) {tetpil=avma;return gerepile(av,tetpil,gcopy(x));}
p3=factor(p4);p1=(GEN)p3[1];p2=(GEN)p3[2];nfact=lg(p1)-1;
av1=avma;p3=cgetg(n+1,17);
for(i=1;i<=nfact;i++)
{
p=(GEN)p1[i];unmodp=gmodulcp(gun,p);fl=1;
while(fl)
{
pk=ker(gmul(unmodp,x));if(lg(pk)==1) fl=0;
else
{
p5=centerlift(pk);p4=gdiv(gmul(x,p5),p);
for(i=1;i<=n;i++) p3[i]=0;
for(j=1;j<lg(p5);j++)
{
j1=n;while(gcmp0(coeff(p5,j1,j))) j1--;
p3[j1]=j;
}
for(j=1;j<=n;j++) if(p3[j]) x[j]=p4[p3[j]];
}
}
if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
}
}
else
{
unmodp=gmodulcp(gun,pp);fl=1;
while(fl)
{
pk=ker(gmul(unmodp,x));if(lg(pk)==1) fl=0;
else
{
p5=centerlift(pk);p4=gdiv(gmul(x,p5),pp);
for(i=1;i<=n;i++) p3[i]=0;
for(j=1;j<lg(p5);j++)
{
j1=n;while(gcmp0(coeff(p5,j1,j))) j1--;
p3[j1]=j;
}
for(j=1;j<=n;j++) if(p3[j]) x[j]=p4[p3[j]];
if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
}
}
}
tetpil=avma;return gerepile(av,tetpil,gcopy(x));
}
GEN matrixqz2(x)
GEN x;
{
long av=avma,tetpil,i,j,k,m,n,fl,lim,in[2];
GEN p1;
if(typ(x)!=19) err(matqzer1);
lim=(avma+bot)>>1;
n=lg(x)-1;if(!n) return gcopy(x);
x=gcopy(x);
m=lg(x[1])-1;
for(i=1;i<=m;i++)
{
do
{
for(fl=0,j=1;(j<=n)&&(fl<2);j++)
if(!gcmp0(coeff(x,i,j))) in[fl++]=j;
if(fl==2)
{
j=(gcmp(gabs(coeff(x,i,in[0])),gabs(coeff(x,i,in[1])))>0)?in[1]:in[0];
p1=(GEN)coeff(x,i,j);
for(k=1;k<=n;k++)
if(k!=j) x[k]=lsub(x[k],gmul(ground(gdiv(coeff(x,i,k),p1)),x[j]));
}
}
while(fl==2);
j=1;while((j<=n)&&gcmp0(coeff(x,i,j))) j++;
if(j<=n) x[j]=lmul(denom(coeff(x,i,j)),x[j]);
if(avma<lim) {tetpil=avma;x=gerepile(av,tetpil,gcopy(x));}
}
tetpil=avma;return gerepile(av,tetpil,hnf(x));
}
GEN matrixqz3(x)
GEN x;
{
long av=avma,av1,tetpil,i,j,j1,k,m,n,fl,lim,in[2];
GEN p1,c,d;
if(typ(x)!=19) err(matqzer1);
lim=(avma+bot)>>1;
n=lg(x)-1;if(!n) return gcopy(x);
x=gcopy(x);m=lg(x[1])-1;
c=cgeti(n+1);for(i=1;i<=n;i++) c[i]=0;
d=cgeti(m+1);for(i=1;i<=m;i++) d[i]=0;
av1=avma;
for(k=1;k<=m;k++)
{
j=1;while((j<=n)&&(c[j]||gcmp0(coeff(x,k,j)))) j++;
if(j<=n)
{
x[j]=ldiv(x[j],coeff(x,k,j));
for(j1=1;j1<=n;j1++) if(j1!=j) x[j1]=lsub(x[j1],gmul(coeff(x,k,j1),x[j]));
c[j]=k;d[k]=j;
}
if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
}
for(i=1;i<=m;i++)
{
do
{
for(fl=0,j=1;(j<=n)&&(fl<2);j++)
if(!gcmp0(coeff(x,i,j))) in[fl++]=j;
if(fl==2)
{
j=(gcmp(gabs(coeff(x,i,in[0])),gabs(coeff(x,i,in[1])))>0)?in[1]:in[0];
p1=(GEN)coeff(x,i,j);
for(k=1;k<=n;k++)
if(k!=j) x[k]=lsub(x[k],gmul(ground(gdiv(coeff(x,i,k),p1)),x[j]));
}
}
while(fl==2);
j=1;while((j<=n)&&gcmp0(coeff(x,i,j))) j++;
if(j<=n) {p1=denom(coeff(x,i,j));if(!gcmp1(p1)) x[j]=lmul(p1,x[j]);}
if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
}
tetpil=avma;return gerepile(av,tetpil,hnf(x));
}
GEN kerint1(x)
GEN x;
{
long av=avma,tetpil;
GEN p1,p2;
p1=matrixqz3(ker(x));
p2=lllint(p1);tetpil=avma;return gerepile(av,tetpil,gmul(p1,p2));
}
GEN kerint2(x)
GEN x;
{
long lx=lg(x), tx=typ(x),i,j,av,av1;
GEN g,p1;
if(tx!=19) err(lller1);
av=avma;
g=cgetg(lx,19);
for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
for(i=1;i<lx;i++)
for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
g=lllgramall(g,1);p1=lllint(g);
av1=avma;return gerepile(av,av1,gmul(g,p1));
}
/*
GEN oldkerint(x)
GEN x;
{
long lx=lg(x), tx=typ(x),i,j,av,av1;
GEN g,p1;
if(tx!=19) err(lller1);
av=avma;
g=cgetg(lx,19);
for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
for(i=1;i<lx;i++)
for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
g=lllgramall0(g,1);p1=lllint(g);
av1=avma;return gerepile(av,av1,gmul(g,p1));
}
*/
GEN kerint(x)
GEN x;
{
long av=avma,av1;
GEN g,p1;
g=lllall0(x,1);p1=lllint(g);
av1=avma;return gerepile(av,av1,gmul(g,p1));
}
GEN intersect(x,y)
GEN x,y;
{
long av=avma,tetpil,i,j,k,p;
GEN z,p1,r;
if((typ(x)!=19)||(typ(y)!=19)) err(interer1);
if((lg(x)==1)||(lg(y)==1)) return cgetg(1,19);
z=ker(concat(x,y));k=lg(x);p=lg(z);
r=cgetg(p,19);
for(j=1;j<p;j++)
{
p1=cgetg(k,18);r[j]=(long)p1;
for(i=1;i<k;i++) p1[i]=coeff(z,i,j);
}
tetpil=avma;return gerepile(av,tetpil,gmul(x,r));
}